#include "kernel/cell.hpp"


#include "yocto/string/vfs-utils.hpp"
#include "yocto/lua/lua-config.hpp"
#include "yocto/math/fcn/zfind.hpp"
#include "yocto/ios/ocstream.hpp"



static Real K1p   = 0;
static Real ratio = 0;
static Real PCO2  = 0;
static AE  *pAE   = NULL;
static NHE *pNHE  = NULL;

static Real L0_org = 0;
static Real Kr_org = 0;
static Real Kt_org = 0;

static Real KAE_org = 0;
static Real gamma_org = 0;

static Real ZpH( Real pH )
{
	const  Real h = Pow(10.0,-pH);
	const  Real b = K1p * (40.0*PCO2/760.0)/h;
	return ratio - pNHE->Sigma(h/pNHE->Kr) / pAE->Sigma(b);
}


static void reset(const Real pH0)
{
	pNHE->L0       = L0_org;
	pNHE->Kr       = Kr_org;
	pNHE->Kt       = Kt_org;
	pNHE->updateC();
	
	pAE->K     = KAE_org;
	pAE->gamma = gamma_org;
	
	PCO2       = 1;
	ratio      = 0;
	ratio      = -ZpH(pH0);
	
}

int main( int argc, char *argv[] )
{
	
	const char *progname = _vfs::get_base_name( argv[0] );
	try
	{
		//======================================================================
		//
		// preparing data
		//
		//======================================================================
		if( argc <= 1 )
			throw exception("Usage: %s config.lua", progname);
		const string config_name = argv[1];
		Lua::State   VM;
		lua_State   *L = VM();
		lua_register(L,"geom_data",geom_data);
		
		Lua::Config::DoFile( L, config_name );
		Cell         cell( L );
		
		//======================================================================
		//
		// initialize steady state
		//
		//======================================================================
		cell.steadyState();
		pAE    = cell.pAE;
		pNHE   = cell.pNHE;
		L0_org = pNHE->L0;
		Kr_org = pNHE->Kr;
		Kt_org = pNHE->Kt;
		
		KAE_org   = pAE->K;
		gamma_org = pAE->gamma;
		
		//======================================================================
		//
		// extract values for computation
		//
		//======================================================================
		const Real K1 = Lua::Config::Get<lua_Number>( L, "K1" );
		const Real KH = Lua::Config::Get<lua_Number>( L, "K_Henry");
		K1p = K1/KH;
		
		
		const size_t  Np    =   64;
		Real          lpmin =  -2.0;
		Real          lpmax =   2.0;
		
		
		zfind<Real>             solve = { 1e-12 };
		numeric<Real>::function zfunc( cfunctor(ZpH));
		const Real pH0 = 7.2;

		//======================================================================
		//
		// new pH as a function of L0/L0,wild_type
		//
		//======================================================================
		
		{
			
			reset(pH0);
			
			
			std::cerr << "Default Ratio=" << ratio << std::endl;
			
			ios::ocstream fp("pH-L0.dat",false);
			fp("#L0/L0_wild pH\n");
			for( size_t i=1; i <= Np; ++i )
			{
				const Real lp = lpmin + ( (i-1) * (lpmax-lpmin))/(Np-1);
				const Real r  = Pow(10,lp);
				const Real L0 = L0_org * r;
				pNHE->L0      = L0;
				const Real pH1 = solve( zfunc,1,15);
				fp("%.7g %.4f\n", r, pH1 );
			}
		}
		
		//======================================================================
		//
		// new pH as a function of Kr,Kt
		//
		//======================================================================
		{
			reset(pH0);
			ios::ocstream fp("pH-KrKt.dat",false);
			fp("#Kr/Kr_wild pH\n");
			for( size_t i=1; i <= Np; ++i )
			{
				const Real lp = lpmin + ( (i-1) * (lpmax-lpmin))/(Np-1);
				const Real r  = Pow(10,lp);
				pNHE->Kr      = Kr_org * r;
				pNHE->Kt      = Kt_org * r;
				pNHE->updateC();
				const Real pH1 = solve( zfunc,1,15);
				fp("%.7g %.4f\n", r, pH1 );
			}
		}
		
		//======================================================================
		//
		// new pH as a function of KAE
		//
		//======================================================================
		{
			reset(pH0);
			ios::ocstream fp("pH-KAE.dat",false);
			fp("#KAE/KAE_wild pH\n");
			for( size_t i=1; i <= Np; ++i )
			{
				const Real lp = lpmin + ( (i-1) * (lpmax-lpmin))/(Np-1);
				const Real r  = Pow(10,lp);
				pAE->K = KAE_org * r;
				const Real pH1 = solve( zfunc,1,15);
				fp("%.7g %.4f\n", r, pH1 );
			}
			
		}
		
		//======================================================================
		//
		// new pH as a function of Hill's coefficient
		//
		//======================================================================
		{
			reset(pH0);
			ios::ocstream fp("pH-alpha.dat",false);
			fp("#alpha/alpha_wild pH\n");
			lpmax = 1.7;
			for( size_t i=1; i <= Np; ++i )
			{
				const Real lp = lpmin + ( (i-1) * (lpmax-lpmin))/(Np-1);
				const Real r  = Pow(10,lp);
				pAE->gamma = gamma_org * r;
				const Real pH1 = solve( zfunc,1,15);
				fp("%.7g %.4f\n", r, pH1 );
			}
			
		}
		
		
		return 0;
	}
	catch( const exception &e )
	{
		std::cerr << e.when() << std::endl;
		std::cerr << e.what() << std::endl;
	}
	catch(...)
	{
		std::cerr << "Unhandled Exception in " << progname << std::endl;
	}
	return -1;
}

#if 0
{
	ios::ocstream fp("pH0.dat",false);
	for( size_t i=1; i <= Np; ++i )
	{
		const Real lp = lpmin + ( (i-1) * (lpmax-lpmin))/(Np-1);
		PCO2 = Pow(10,lp);
		std::cerr << "PCO2/40mmHg=" << PCO2 << std::endl;
		for( size_t j=1; j<=Nr; ++j )
		{
			const Real lr = lrmin + ( (j-1) * (lrmax-lrmin) )/(Nr-1);
			ratio = Pow(10,lr);
			std::cerr << "\tV_AE/V_NHE=" << ratio << std::endl;
			Real pH_lo = 1;
			Real pH_hi = 15;
			const Real pH0 = solve( zfunc, pH_lo, pH_hi);
			std::cerr << "\t\tpH=" << pH0 << std::endl;
			fp("%.4f",pH0);
			if( j < Nr ) fp(" ");
		}
		fp("\n");
		}
		}
#endif
		
		
